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I. INTRODUCTION 

Two areas at the forefront of research in Bose-Einstein 
condensates (BECs) over the last few years have been 
optical lattices and the hyperfine degree of freedom [1, 
2]. Optical lattices allow one to explore, in both the 
mean- field and quantum regimes, the effects of periodic 
potentials on bosons. These studies complement the vast 
body of knowledge concerning the behavior of fcrmions 
in periodic potentials coming from solid state physics. 
Optical lattices are free of defects and disorder and the 
atomic potential has a simple closed form. In contrast, 
an electron in a solid is subject to a complex, imperfectly 
known potential that is usually marred by defects. 

The hyperfine states of the atoms making up a BEG 
allow one to construct exotic spin structures based on 
the occupation of different hyperfine spin states of the 
form \F, mp)- BECs of this kind, which are called spinor 
or multicomponent, have a vector order parameter. Re- 
search on multicomponent BECs has been instrumental 
in reaching important milestones in BEC research, in- 
cluding the creation of a quantum vortex and the sub- 
sequent demonstration that a BEC made from a weakly 
interacting alkali gas is supcrfluid [3, 4]. Recently, exper- 
imentalists have placed multicomponent BECs in optical 
lattices [5, 6]. 

In this article, we construct exact solutions to the 
mean-field equations of motion for multicomponent 
BECs in one, two and three dimensional optical lattices. 
Band theory was invented as a tool to analyze station- 
ary solutions of the linear Schrodinger equation with a 
periodic potential, a problem that does not have a so- 
lution in closed analytic form for realistic potentials. In 
contrast, for the nonlinear Schrodinger equation (NLS) 
with a sinusoidal optical potential, exact, closed- form 
stationary solutions have been discovered [7-14]. Here 
we generalize and extend previous work to an overar- 



ching and rigorous treatment of certain classes of exact 
solutions describing the dynamics of condensates with s 
components in D dimensions. Our formalism permits us 
to construct exact, nonstationary solutions of the vec- 
tor NLS. This article brings together our new work and 
past treatments into a single, general, rigorous analyti- 
cal framework. At the same time, we elucidate the most 
experimentally relevant and aesthetically pleasing solu- 
tions. We also perform a full nonlinear stability analysis 
where computationally tractable. A number of surprising 
stability regimes present themselves, as we will demon- 
strate. 

The basic idea that leads to our exact solutions is to 
use the nonlinearity to cancel the spatial variation of the 
potential, leading to an effective free particle problem. 
(Clearly, this is not possible for the linear Schrodinger 
equation.) This idea is due to Bronski, Carr, Deconinck, 
and Kutz, who originally applied it to a one-component, 
or scalar, BEC in a one-dimensional periodic potential [7- 
9]. Deconinck, Frigyik and Kutz extended this work on 
one-component BECs to higher dimensions [10, 11], and 
to multicomponent condensates in one dimension [13], 
but not to multicomponent BECs in two and three di- 
mensions. Although the higher-dimensional Jacobi el- 
liptic periodic potentials considered by Deconinck et al. 
do not include the square, rectangular and simple cubic 
optical potentials readily available in the lab, Hai and 
coworkers have shown that the cancellation technique can 
be used to construct solutions for a scalar condensate on 
a square optical lattice [12]. 

Our extension of the cancellation technique to conden- 
sates with an arbitrary number of components s in a sin- 
uoidal optical potential of arbitrary dimension D leads to 
an enormous number of new solutions, many of great ex- 
perimental import. A crucial and challenging step in the 
cancellation technique is to find a solution ansatz that, 
for a given potential, allows the cancellation to occur. At 
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the same time, the solution ansatz must satisfy the free 
particle Schrodinger equation. A key aspect of our work 
is the introduction of a novel, very general solution ansatz 
that allows the solution of a wide range of problems. 

Mean-field theory, which for a single-component con- 
densate takes the form of the scalar NLS or Gross- 
Pitaevskii equation [15], has been quite successful in de- 
scribing experiments on multicomponent BECs in optical 
lattices [5, 16]. However, there are important approxima- 
tions underlying its use. First, the tunneling or hopping 
energy th must be much larger than the on-site interac- 
tion energy U so that the system is far from the Mott 
insulating regime [17, 18]. This means that the poten- 
tial barriers are not so high that the sites lose mutual 
phase coherence, and that a full many body quantum 
Fock state treatment is not needed [20]. Second, three- 
body and other loss processes are neglected. Third, quan- 
tum fluctuations are ignored, as is always the case when 
the NLS is applied to BECs [15, 20]. Fourth, any possible 
resonances induced by the lattice or dimensional confine- 
ment arc neglected [21]. Fifth, when treating D = 1 
and D = 2, mean- field theory requires that the confining 
potential that reduces the effective dimensionality have 
a length scale smaller than or of the order of the heal- 
ing length but larger than the scattering length [22-24], 
i.e., the underlying scattering process must remain three 
dimensional. 

Three important conditions must be met if our can- 
cellation technique is to apply. First, we cannot include 
the low frequency harmonic trap often, but not always, 
present in experiments. Such a trap is used to keep atoms 
from spilling off the edge of a finite lattice. We require 
the potential to be sinusoidal, although we do allow the 
lattice constant to be different in each direction, leading 
to, e.g., a rectangular lattice in two dimensions. Sec- 
ond, for condensates with three or more components, the 
mean-field theory normally includes coherent couplings 
between different components of the vector order param- 
eter; for two components, such couplings are prevented 
by angular momentum selection rules [25, 26]. We require 
incoherent couplings only, which is appropriate for s = 2. 
For s > 2 our treatment is always correct for sufficiently 
short time scales [27]. Our treatment can also be cor- 
rect at arbitrary times when the hyperfine components 
arc chosen from separate manifolds F. Finally, we can- 
not treat a mixture of scalar BECs of different masses, 
despite the fact that such a system can be described by a 
vector mean-field theory with incoherent couplings only. 

Although certain conditions must be satisfied for it to 
be applicable, our cancellation technique yields a panoply 
of exact solutions of great physical significance. For ex- 
ample, for a two-component condensate on a rectangular 
optical lattice, we find temporally periodic solutions in 
which the optical lattice is divided into two sublattices, 
and the condensate components oscillate back and forth 
between these sublattices. For the square optical lattice, 
we find a vortex-anti-vortex array for a scalar conden- 
sate, while for two-component condensates we obtain ex- 



otic solutions in which the optical lattice is divided into a 
total of four sublattices, and the condensate components 
move cyclically between these sublattices. As the dimen- 
sion D and the number of components s are increased, 
the number of solutions our technique generates grows 
rapidly. The number of solution types is so vast in three 
dimensions (3D) that, for the sake of brevity, we limit 
our discussion to stationary solutions with a high degree 
of symmetry and to two examples of non-stationary so- 
lutions. 

The article is organized as follows. In Sec. II, we intro- 
duce the mean-field equations of motion and the optical 
potentials we will study. In Sec. Ill, a complete and rigor- 
ous treatment of exact dynamical solutions for arbitrary 
dimensions D and number of components s is presented. 
In Sec. IV, we treat select cases in detail, giving exam- 
ples of how to apply the results of Sec. Ill; of particular 
interest are the vortex-anti-vortex array we find in 2D 
and the array of intersecting vortex lines we find in 3D. 
In Sec. V, we present detailed stability studies of im- 
portant solution classes, including the vortex-anti- vortex 
array, and an explicit connection to experimental units. 
Finally, in Sec. VI, we conclude. 



II. THE MEAN-FIELD EQUATIONS OF 
MOTION 

Consider an s-component BEC in D dimensions with 
incoherent couplings between components. The conden- 
sate is subject to an optical potential formed by D lin- 
early polarized retrofiected light waves. Let k; be the 
wave vector of the l^^ light wave, where I G {1, . . . , D} 
will be called the directional index. We will restrict our 
attention to optical lattices with k/ • k;/ =0 for / 7^ 
In particular, we will study BECs in one-dimensional, 
square, rectangular, and cubic optical lattices. For a 
review of optical potentials for neutral atoms, see, for 
instance, Ref. [28]. 

The mean-field equations of motion are 



ih 



dt 



2m 



^J, (1) 



where ^j(r, t) is the j component of the conden- 
sate order parameter, j g {l,...,s}, and the posi- 
tion r is a vector in D dimensions. The j^^ compo- 
nent of the order parameter may be written ipjir^t) = 
^ rij (r, t) expliSj (r, t)], where nj(r, t) is the number den- 
sity of the j**^ component at position r and Vj(r,t) = 
{h/m)'V Sj{r,t) is its velocity at that point [15]. Atoms 
in the j*'^ component are subject to the optical potential 
Vj . The coefficients gjji of the nonlinear terms describe 
the binary interaction of an atom in component j and 
an atom in component j': explicitly, gjji = 47r?i^ ajj'/m, 
where ajji is the s-wave scattering length and m is the 
atomic mass, which, as stated in Sec. I, is assumed to 
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be independent of the component index j. Note that the 
nonlinear coefhcient gjji is renormahzed by transverse 
confinement as briefly alluded to in Sec. I for D = 1 and 
D = 2; see the references given there for more details. 
The nx n real, symmetric matrix M = {gjj'} will be re- 
ferred to as the interaction matrix. We will assume that 
all of the diagonal elements of M are nonzero, as is the 
case in experiments. 

The atoms of the j**^ component are subject to the 
optical potential 

1 ^ 

T/,(r) = --p,^e2cos2(k,.r), (2) 
;=i 

where pj is the atomic polarizability of an atom in the j*'' 
hypcrfinc state and ei is the electric field amplitude of the 
l^^ standing light wave. The wave vector k; determines 
the lattice constant in the l^^ direction. For convenience, 
we let Vji = jPjcf. Equation (2) then becomes 

D 

V,{r) = -J2v,iC0s\ki-r). (3) 

1=1 

We assume that all of the Pj's and ej's are nonzero, so 
that none of the Vji^s vanish; if Vji were zero, the atoms 
of the j*"^ component would be subject to an optical po- 
tential independent of the l^^ spatial coordinate. Note 
that the atomic polarizability pj generally depends on 
the component, or hypcrfine state, j. 

To illustrate how one arrives at the potential given by 
Eq. (2), consider the case D = 2. The total electric field 
E = El -I- E2, where 

E/ = e/ cos(k/ • r + X;) cos(c|k; |i + 0/) . (4) 

The vector e; has constant, real components and is or- 
thogonal to k/ , and c is the speed of light in vacuum. By 
shifting the location of the origin if necessary, one can 
arrange for both of the spatial phases xi to vanish. The 
temporally averaged intensity is then 

/ = ie?cos2(ki •r)-Hie^cos^(k2 •O-^Jfc.^fc, 

xei • 62 cos((/)i — cos(ki • r) cos(k2 • r) . (5) 

The j*'^ component of the BEC is subject to the external 
optical potential Vj = —\pjl. The optical potentials Vj 
take the form (2) with _D = 2 if and only if the term 
in / coming from the interference of the two light waves 
vanishes. If ki and k2 differ, the interference term is zero 
and the optical potential has rectangular symmetry. If 
fci = ^2, on the other hand, the interference term van- 
ishes if either ei • 62 = or cos(0i — 4>2) = 0. The optical 
lattice is then simply a square lattice. Finally, if fci = k2 
and < |ei • 62 cos(0i — ^2)! < 6162, the structure of the 
optical lattice is more complex [29] . 

The interference terms in the intensity can also be 
made to vanish for D = 3: for example, we can choose 



the vectors to be orthogonal to one another. If 
ki ~ k2 = k^, the optical lattice has simple cubic sym- 
metry. 

For simplicity, in this paper we will confine ourselves to 
optical potentials in which the interference terms vanish. 
The potential is then given by Eq. (2). It is worth noting, 
though, that our solution techniques can be generalized 
to optical potentials with nonzero interference terms. 

III. METHODS OF CONSTRUCTING EXACT 
SOLUTIONS TO THE MEAN-FIELD 
EQUATIONS OF MOTION 

The time evolution of the order parameter is described 
by the mean-field equations of motion (1) with the optical 
potentials (3). We will seek solutions to this problem in 
which each of the effective potentials 

s 

U,{r,t)^V,{r) + Y,9jj'\i^j'ir,t)\' (6) 
i'=i 

is constant. Solutions of this kind will be referred to as 
potential- canceling (PC) solutions and our solution tech- 
nique will be called the cancellation method because for 
each j € {l,2,...,s}, the spatial variation of the op- 
tical potential Vj{r) is canceled by the variation of the 
term X]/=i .9j j' IV'j' (r, t)P, rendering the effective poten- 
tial Uj constant. PC solutions reduce the coupled non- 
linear mean- field equations of motion (1) to uncoupled 
linear Schrodinger equations with constant potential: 

for j€{l,2,...s}. 

For the JTj's to be constant, the tpj^s must be linear 
combinations of terms that vary sinusoidally with posi- 
tion. To be precise, we seek solutions of the form 

D 

V,(r,t) = e-'"^*^^,zcos(kz • r)e-'"'*, (8) 
/=o 

where huji — Ji^kf /2m is the energy of a free particle of 
mass m with wave number ki in the absence of nonlin- 
carity. Note that the possible values of the directional 
index I have been extended to include ^ = in order 
to simplify the notation: ko = so that the I = term 
gives rise to a constant offset in the order parameter. The 
coefficients Aji in our solution ansatz (8) are in general 
complex, while the frequencies Qj are real. The Aji^s are 
constrained by the requirement that the effective poten- 
tials Ui, U2, ■ ■ ■ ,Us are constant. These constraints will 
be discussed in detail below. 

Substituting Eq. (8) into Eq. (1), we see that 

s 

Vj = -J2an'\^r\^ + ^^j, (9) 
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for J £ {1, 2, ... s}. This means that the effective poten- 
tial Uj takes on the constant value hflj. Inserting Eq. (8) 
into Eq. (9) and comparing the resulting expression for 
Vj with Eq. (3), we obtain 

D D D / s \ 

^ cos2(k, • r) = ^ ^ gn'AriA*,, 

1=1 1=0 l'=0 \j' = l / 

X cos(kz • r) cos(kr • r)e-*'"'-"''^* - hn^, (10) 

for j £ {1,2,..., s}. Equation (10) yields a set of alge- 
braic equations that the coefficients Aji must satisfy. The 
first set of equations ensure that the cross terms on the 
right hand side of Eq. (10) vanish. Specifically, for each 
pair of integers {I, I') with < I < I' < D, one obtains a 
set of conditions. If uji ^ uji', we must have 

s 

Y,9n'AriA*,i,^0 (11) 

for j G {1, . . . , s}. On the other hand, if loi = loii , it is 
sufficient to impose the weaker conditions 

^[Y.gn'A,nA*,)j=Q. (12) 

Equating the coefficients of the terms that are propor- 
tional to cos^(k; • r) on either side of Eq. (10), we see 
that 

s 

Vji^Y.9n'\Ari? (13) 

for j e {1, 2, . . . , s} and ? e {1, 2, . . . , D}. (Note that this 
equation does not apply for / = 0.) Finally, the constant 
term on the right hand side of Eq. (10) must vanish, and 

so 

for j e {l,2,...,s}. 

The task of finding solutions to the coupled nonlinear 
partial differential equations (1) has now been reduced 
to solving a system of algebraic equations: Eq. (13) must 
be solved for the coefficients Aji subject to the condi- 
tions (11) or (12) for each pair (/,/') with I < I'. A 
solution to these equations does not necessarily exist. If 
a solution does exist, Eq. (14) yields the frequencies flj. 

The solution, if it exists, is not uniquely specified by 
the system of algebraic equations. To see this, let 

Ai = {Au,...,A,if (15) 

for / e {0, 1, 2, . . . , £)}. Equation (13) determines the 
norm of the vector A; for I G {1,2,. ..,£)} but does not 
constrain the magnitude of Aq. In addition, Eqs. (11) 



and (12) with 1 = place constraints on the direction of 
Aq but not its norm. The quantity |Aop is therefore a 
free parameter. 

The length of the vector Aq is determined if the spa- 
tial average of the total density is given, as we will now 
establish. At time t, the number density of the j^^ com- 
ponent at position r is nj{r,t) = |-0j(r, t)p. Let (/) 
denote the spatial average of an arbitrary function /(r). 
Using Eq. (8), we find that the spatial average of the 
total number density 

s 

(n) ^ ^(n,) (16) 

is given by 

(n) = |Aop + if]|A,p. (17) 

1=1 

Equation (17) has a solution for |Aop if and only if 

D 

2{n)>Y,\A,f; (18) 
1=1 

if this condition is met, |Aop is uniquely specified. 

Equation (11) states that the vector 
{AiiAli,,...,AniA^i,)'^ is in the kernel of M, while 
if Eq. (12) applies, the real part of this vector must be 
in the kernel of M. For this reason, most (but not all) 
of the solutions we obtain will be for cases in which the 
atomic interactions are such that det M = 0. 

We will now consider three particularly interesting and 
physically important special cases in which the algebraic 
conditions that must be solved to yield a solution sim- 
plify dramatically. These special cases will be referred 
to as Special Cases A, B and C. We will also provide an 
example that shows that the formalism just developed 
yields solutions to the mean-field equations of motion 
even when the interaction matrix is nonsingular. This 
example appears in Subsection HID. 

A. Factorizable Equations of Motion 

A particularly simple special case is obtained when the 
rank of M is unity. This is true to an excellent approxi- 
mation for the two-component condensates first produced 
by the JILA group that consist of two different hypcrfine 
spin states of ^^Rb: gn, gi2 and 522 arc known to the 
1% level, and arc in the proportion 1.03 : 1 : 0.97 [30]. 
As a result, detM/TrM is zero to within experimental 
error. 

Because M has rank 1 and is a symmetric matrix, there 
are nonzero, dimcnsionlcss, real numbers Aj and a cr = 
±1 such that 

933' = crgAjAjv (19) 
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for all j and j' . The quantity g is a positive constant with 
dimensions of energy times volume, which is inserted in 
Eq. (19) to render the A^ 's dimensionless. The magnitude 
of g is arbitrary but fixed and, if desired, may be taken 
to be the typical magnitude of the interaction coefficients 

9jj' ■ 

Let L = (Ai, A2, . . . , As)"^. Equation (19) may then be 
written 

M^agLlT^, (20) 

showing that when the rank of the interaction matrix is 
1, M can be factored. 

Let A be the s x s matrix with elements 

A,,' ^ \j5,,r ■ (21) 

For each pair of directional indices {I, I') with < I < 
I' < D, Eq. (11) reduces to the single condition 

A|AAj, =0, (22) 

while Eq. (12) becomes 

3fJ(A|AAzO = . (23) 

The relations (13) are now 

^p,ef = agX.AjAAi , (24) 

where j G {1, . . . , s} and / G {1, 2, . . . , D}. Equation (24) 
has a solution if and only if 

Pj^pXj, (25) 

for all J, where p is a nonzero real constant. If this is the 
case, then 

Vj = X,V, (26) 

where 

1 1 " 

V ^ --pi ^--pY^ef cos' {ki-v). (27) 

1=1 

Equation (24) then becomes 

agAjAAi = -^pe} ; (28) 

this holds for Z G 1, 2, . . . , _D. Equation (14) shows that 

= Ajl7, (29) 

where 

f] = ^aJAAq. (30) 

Finally, Eq. (9) reduces to 

V ^-ogil^^Aii + m, (31) 



where 

^ = (^-1,...,^.)^ (32) 

is the vector order parameter. 

As before, |Aop is a free parameter unless an addi- 
tional constraint is applied. If the spatial average of the 
total density (n) is given and 

(n)>lf]|A,p, (33) 

1=1 

then 

|AoP = (n)-^E|A,p. (34) 

i=\ 

If M has rank 1 and the condition (25) holds, we call 
the mean-field equations of motion factorizable. For this 
case, which we will refer to as Special Case A, the equa- 
tions of motion (1) assume the simpler form 

h' 

in-^= V V + crg(i/'Ui/>)A^/> + yAi/> . (35) 

ot 2m 

To find exact solutions to the coupled nonlinear partial 
differential equations (35), Eqs. (28) and (34) must be 
solved for Aq, Ai, . . . , A^i subject to the conditions (22) 
or (23) for each pair {1,1') with < I < I' < D. If 
a solution to these equations exists, Eqs. (29) and (30) 
yield the frequencies flj. 

In Appendix A we demonstrate that if the mean-field 
equations of motion are factorizable, an additional sim- 
plification can be made: without loss of generality, all of 
the Aj 's may be taken to be of unit modulus. Therefore, 
for the remainder of the paper, when we discuss Special 
Case A, we will assume that each of the Aj's is equal to 
±1. 



1. Constructing Exact Non-stationary Solutions Using a 
Transformation 

If we have a PC solution ^0 to the equation of mo- 
tion (35), under certain circumstances we can construct 
new solutions by transforming xp. Let P be an invertible 
s X s matrix, and suppose that is given by Eq. (8) and 
satisfies Eq. (35). We set 

C(r,i) =T(P)V'(r,0, (36) 

where 

T(P) = exp{-inAt)Pexp{inAt). (37) 

From Eq. (8), it follows that if 

pUp = A, (38) 

then C(r,i) is also a solution to the equation of mo- 
tion (35). An important aspect of this transformation 
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is that even if ■0(r, t) is stationary, ^(r, t) can turn out to 
be nonstationary. Thus, by transforming a single solu- 
tion ■0, we obtain a set of stationary and nonstationary 
solutions. We will call each such set a P-set. 

For a given A, let >S'(A) be the set of invertible s x s 
matrices P that satisfy Eq. (38). It is straightforward 
to show that the set of matrices T(P) with P G 5 (A) 
forms a group. We will call this group by G{A). Later 
in the paper we will study two examples in which G{A) 
is a continuous group, i.e., it has an uncountably infinite 
number of elements. As a result, the P-set is uncountably 
infinite in these examples. 

B. Factorizable Equations of Motion with Equal 
Atomic Polarizabilities 

A particularly important factorizable problem has 
Xj — 1 for all J, so that A is the identity matrix X. In this 
case, which we will call Special Case B, all of the inter- 
action strengths Qjj' have the value ag, and the atomic 
polarizabilities pj are all equal. 

The equation of motion (35) with A = T and V ~ 
was first studied by Manakov [31] and is now known as 
the Manakov equation. We will extend this terminology 
by also calling Eq. (35) with A = X and nonzero potential 
V the Manakov equation. 

The Manakov Case, i.e.. Special Case B, is of consid- 
erable physical interest. Provided that the atoms are not 
too close to resonance, the pj's are to a good approxi- 
mation equal [32]. The interaction strengths are nearly 
equal in two-component ^''Rb condensates [30, 33]. As a 
result, the dynamics of these condensates are reasonably 
well described by the Manakov equation with s = 2. 

Three-component ^'^Na condensates with hyperfine 
spin P = 1 were first studied by the MIT group [34, 35]. 
For F = 1 spinor condensates, the interaction strengths 
Qjji are identical and there are no incoherent couplings 
if Iq and I2 are equal, where lyr is the s-wave scattering 
length for two colliding atoms with total hyperfine spin 

[25, 26]. Since the difference I2 — Iq is small compared 
to Iq for ^^Na [35, 36], it is a reasonable approximation 
use the Manakov equation with s = 3 to model the three- 
component condensates produced by the MIT group, at 
least for the initial stage of the time evolution. 

For the Manakov case, Eq. (28) reduces to 

\M'-^el (39) 

where / £ {1,2, . . . , D}. This shows that ap must be 
positive for there to be a solution to Eq. (39). Thus, for 
the remainder of the paper, when we discuss Special Case 
B, we will assume that ap > 0. For convenience, let 

for I = 1,2, . . . , D. Equation (39) is then simply ]A;] = 
ai. 



The optical potentials Vj are all equal to V for Special 
Case B. Equation (31) shows that 

V = -agn + nn, (41) 

where n = is the total condensate number density. 
The total density is independent of time and varies si- 
nusoidally with position. The maxima of n are located 
at the potential minima for cr = -1-1. In contrast, for 
CT = — 1, the maxima of n are located at the maxima 
of the potential. This leads to an obvious instability, as 
pointed out by Bronski et at [9] for the single-component 
case. Accordingly, for the remainder of the paper, we will 
limit our attention to the case a = +1 whenever we study 
a Case B problem [37]. Since we have already assumed 
that ap > 0, this means that for Case B the atomic po- 
larizability p will be taken to be positive throughout the 
remainder of the paper. 

The condition (38) is particularly simple for the Man- 
akov equation: P can be any unitary matrix. Since A is 
the identity matrix, T{P) = P and C is a unitary trans- 
formation of xp. Unitary transformations of solutions to 
the Manakov equation with an external potential have 
been studied elsewhere in one spatial dimension [14, 38]. 
The transformation given by Eqs. (36)~(38) generalizes 
that work to problems in which the Xj 's are not all iden- 
tical, as well as to higher spatial dimensions. 

C. Factorizable Equations of Motion for 
Two-Component Condensates with pi = —p2 

Consider a two-component BEC with factorizable 
equations of motion. Recall that Ai and A2 have unit 
modulus and are real. As a result, there are four pos- 
sibilities: (i) Ai = A2 = 1, (ii) Ai = A2 = —1; (iii) 
Ai = — A2 — 1 and (iv) — Ai A2 = 1. Case (i) has al- 
ready been discussed: it is the Manakov case with s = 2. 
The equations of motion for case (ii) are unchanged if 
we reverse the signs of Ai, A2 and p, and so case (n) is 
identical to case (i). In precisely the same way, case (iv) 
is equivalent to case (iii). In this section, we will study 
case (iii). 

The case in which a two-component BEC with factor- 
izable equations of motion has Ai = — A2 ~ 1 will be 
referred to as Special Case C or the factorizable-with- 
opposite-polarizabilities (FOP) case. In this case, the 
atomic polarizabilities pi and p2 have opposite signs and 
the interaction matrix 

M^ag(^^^ "J). (42) 

If a is positive, atoms in the same condensate compo- 
nent repel each other and atoms in different condensate 
components attract. The situation is reversed if a is neg- 
ative. Finally, note that by switching the labels of the 
two components if necessary, we can arrange for ap to be 
negative. We will always take ap to be negative when we 
discuss Special Case C. 
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For Special Case C, the invertible matrix 

A N ( cosh A — sinhA\ 
^(^)^(-sinhA cosh A j' (43) 

satisfies the condition (38) for arbitrary real A. Equa- 
tions (36) and (37) with P = P(A) and real A therefore 
defines a P-set. 



D. Two-Component Condensates with a 
Non-singular Interaction Matrix 

In the special cases discussed so far, the interaction 
matrix M was singular. The cancellation method, how- 
ever, does yield solutions even if det M is nonzero. To 
illustrate this point, let us consider a two-component con- 
densate in one dimension with det Af 7^ 0. 

It follows from Eq. (11) that 



(44) 



A\\ and Ai\ cannot both vanish because the potential 
coefficients Vj\ are nonzero. If both A\\ and A21 are 
nonzero, Aio = A20 = 0, f2i = ^2 = 0, and 



V'j = Ajx cos(A:ia;)e~"^i* 



(45) 



for j = 1,2. A solution of this form exists only if the 
equations 



(46) 



have a solution for A\\ and A^i. The equations (46) have 
a solution if and only if 



and 



XI = (522^11 - 512^^21)/ det U > (47) 



X2 = (311^^21 - .921^^11)/ det M > 0. (48) 



If the inequalities (47) and (48) hold, Aji = y^^. yields 
a solution. 

We next turn to the case in which only one of the 
Aji^s is nonzero. The equations (46) have a solution with 
All = if 



Pi ^ P2 
912 922 



> 0. 



(49) 



On the other hand, Eqs. (46) have a solution with A21 
if 



311 521 



(50) 



If the condition (50) is satisfied, we can switch the label- 
ing of the two condensate components, yielding Eq. (49). 
It is therefore sufficient to consider the case in which the 



condition (49) holds. In this case, A21 = \/ 1^11/512 and 
A20 = All = 0. It follows that tpi = Aioe~'"i* and 
-02 = ^216"*'^"^+"^'* cos(fcia;). Equation (14) becomes 
hflj = gjijAiop, and, without loss of generality, we may 
take AiQ to be real. We conclude that the two-component 
order parameter is given by 



aoe 



(51) 



and 



?/>2 = J^e-'("i+f=^i°"/'''*cosfcia;, (52) 

V .912 

where ag = Aiq is an arbitrary real constant. 



IV. APPLICATION OF ANALYTICAL 
TECHNIQUES TO SELECT CASES 

We will now construct solutions to the mean-field equa- 
tions of motion (1) using the exact analytical methods 
developed in the preceding section. We will start with 
the simplest cases as an introduction to the application 
of our solution methods, and to make connections with 
the relatively simple solutions to be found in the litera- 
ture. We will then move on to progressively more rich 
and complex problems with higher dimensions and/or 
more condensate components than have previously been 
considered. 



A. Solutions on a One Dimensional Optical Lattice 

It is convenient to orient the x axis along ki, so that 
ki ~ kix. Since ujq ^ toi, Eq. (11) applies with / — 
and /' 1. 



1. One-Component Condensates 

We will begin with the simplest case, D = s = 1. For 
s = 1, Eq. (11) reduces to giiAii^A\i = 0. It follows 
that Aio and/or An must vanish. Equation (13) reduces 
to Vii ~ 311 1 All p. Since Vn has been assumed to be 
nonzero. An cannot vanish, and hence Aiq = 0. Equa- 
tion (14) then shows that VLi = 0. There is a solution of 
the form of (8) only if gn and Vn have the same sign. If 
this is the case, we have the solution 



^1 




(53) 



previously found by Bronski et al. [8, 9]. Note that there 
is a PC solution only if the spatially averaged condensate 
density {rii) happens to be Vii/(25ii). Although this is a 
very restrictive condition, this simple case is nevertheless 
useful in the development of nonlinear band theory [39] . 
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2. Two- Component Condensates 

We only found a single stationary solution for a one- 
component condensate in ID. For two-component con- 
densates of the Manakov and FOP types, we find a much 
larger parameter space of solutions, including nonstation- 
ary solutions. 

We briefly touched on solutions for two-component 
condensates in one dimension in Section HID. In ob- 
taining the solution given by Eqs. (51) and (52), we did 
not use our assumption that dct M ^ 0. Moreover, the 
condition (49) holds for both case B and case C. As a 
result, the solution is also valid for cases B and C. In 
both cases, the solution takes the form = where 

^(1) ^ ^-^.galAt/n f «0 \ (54) 

cos(Kia;)e ^ y ^ ' 

The solutions for s = 2 constructed to this point are 
stationary, and have previously been obtained by Decon- 
inck et al. [13]. Let us now consider the Manakov and 
FOP cases B and C. We will demonstrate that for these 
two cases there are nonstationary solutions in the same 
P-set as Eq. (54), where a P-set is the set of solutions 
connected by a matrix transformation of the order pa- 
rameter as described in Sec. Ill A 1; in the Manakov case, 
the transformation is just a unitary transformation. 

For Case B, Eq. (22) becomes Aq-AJ = 0. By changing 
the phase of ipj if necessary, we can arrange for Ajo to 
be real for j = 1,2. We can arrange for An to be real by 
changing the zero of time if needed. It then follows that 
A21 is real as well, and so the vectors Aq and Ai have real 
components. Recalling that |Ai| = ai, we obtain Ag = 
ao(cos6', sin6')-^ and Ai = ai(— sin^, cos 0)^, where qq 
and 9 are arbitrary real constants. The corresponding 
order parameter is 

i/'(r,t) =e-*f"o*/'"'[Ao + AiCos(fcia;)e-''"i*] (55) 

by Eqs. (8) and (30). Recasting this solution, we have 
i/'(r,t) ^ P{0)ipi^\r,t), where 

p{e)^h'^, -''""l) (56) 

^ ' \^sm^ cosOJ ^ ' 

is a unitary matrix. The solution i/'(r, t), which has been 
previously described by Bradley et al. [14], is therefore a 
unitary transformation of the stationary solution (54). If 
sin(20) is nonzero, it is a nonstationary solution because 
the condensate component densities 

ni = ^ alecs'^ 9 + alsin^ 9 cos'^{kix) 

—aoai sm{29) cos(fcia;) cos(wit) (57) 

and 

n2 = |V'2|^ = flo^in^ 6* a^cos^ 0cos^(/i;ia;) 

+aoai sin(20) cos(fcia;) cos(wit) (58) 



oscillate in time with period T = 27r/a;i. We conclude 
that by performing unitary transformations of the single 
solution tl^i^^ , we generate an uncountably infinite P-set 
of stationary and nonstationary solutions. 

What is the physical meaning of the solution (55)? The 
external potentials Vi and V2 coincide and are equal to 
V = ~\pe\ cos^(fcia;), and so the potential minima occur 
at the points x = qir/ki, where q is any integer. We di- 
vide the lattice of potential minima into two sublattices: 
sublattice 1 with even q, and sublattice 2 with odd q. 
The total condensate density n = ni + n2 does not vary 
in time, and its maxima occur at the points x = qn/ki, 
where q is any integer. Suppose for the sake of speci- 
ficity that ao is positive and that n /2 < 9 < n. At time 
< = 0, the maxima of ni are on sublattice 1, while at time 
t = T/2, the maxima of ni are on sublattice 2. At time 
t = T, the maxima of ni are again on sublattice 1. The 
maxima of n2 also oscillate between sublattices 1 and 2, 
but the oscillations of 712 lag those of ni by half a period, 
ensuring that n is time- independent. 

The spatial average of the total number density (n) 
is Gq + ^af. If (71) is given, there is a solution of the 
form (55) if and only if (n) > ^a^. In constrast to the 
single-component case, we obtain a solution not just for 
a single value of (n) , but for a whole rangle of (n) values. 

This concludes our discussion of the relation between 
the solutions obtained using the general formalism of 
Sec. Ill and the existing literature for one dimension. 
Let us now turn to the novel FOP case. Special Case C. 
We can again arrange for the vectors Aq and Ai to be 
real. Equations (22) and (28) have the solution Aq = 
ao(cosh A, — sinh A)'^ and Ai = ai(— sinh A, cosh A)-^ 
valid for arbitrary real ao and A. Since hfl = aga^, the 
condensate component order parameters are 

= e"""S"«*/'' [aoCOshA 

-aisinhAcos(fcia;)e-"i*] (59) 

and 

^2 = e^'^sa^*/^ [-aosinhA 

-haisinhAcos(fcia;)e""i*] . (60) 

In contrast to the solution (55) we constructed for the 
Manakov case, the component densities ni and 712 oscil- 
late in phase between the two sublattices, and the total 
density n oscillates in time as well. Since the spatial aver- 
age of the total number density (n) = cosh(2A)(ao-|-iaf ), 
we obtain solutions provided that (n) >\a\. 
The vector order parameter may be written 

= cxp{-inAt)P{A) exp{inAt)ipi^^ , (61) 

where P(A) is defined by Eq. (43). Once again, the non- 
stationary solution can be constructed by transforming 
the stationary solution (54) and we have an uncountably 
infinite P-set. 
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3. Three-Component Condensates 

We will not carry out the analysis for all possible cases 
for three components. However, we will touch on two new 
features that arise when we go from two components to 
three. 

For two-component condensates governed by the Man- 
akov equations of motion, we showed that the cancella- 
tion method only yields solutions in which the oscillations 
of the two components between the sublattices are 180° 
out of phase. Three-component condensates have an ad- 
ditional degree of freedom, and this leads to solutions 
with a wide range of relative phases. 

We will restrict our attention to Case B and to so- 
lutions with AiQ = A20 = A^Q > and nonzero coef- 
ficients Aji. By changing the zero of time if needed, 
we can arrange for An to be real and positive. Set 
Aji = |^ji|e* for j = 1,2,3 and note that (/>! = 0. 
The condition (22) gives An + A21 + A31 = 0. Recall 
that I Alp = An+A^i+Ali ~ af. Clearly, there are so- 
lutions in which both A21 and A31 are real. In solutions 
of this type, two of the components oscillate in phase 
with one another and the other component is 180° out of 
phase. 

In addition to these solutions, there is a solution for 
any 02 and ^3 satisfying the conditions 

< (/)2 < TT < 03 < 27r (62) 

and 

< 03 - 02 < TT . (63) 

The condensate component densities are 

|V,f = A% + \A,,\'cos\hx) 

+2Ajo\Aji\ cos(fcix) cos(ti;i< — (pj). (64) 

Thus, a wide variety of phase relationships among the 
condensate component densities are possible for s = 3. 
If s is increased still further, the range of possible phase 
relationships grows rapidly. 

For three-component condensates, it is possible for the 
interaction matrix M to be singular and to have rank 
greater than one. Even though the equations of motion 
are not factorizable, the cancellation method developed 
at the outset of Section III can be applied to yield so- 
lutions for problems of this type. We will illustrate this 
with a one-dimensional example. 

Suppose the eigenvalues of M are Ti 7^ 0, T2 7^ and 
T3 = 0, and let the fii's be the associated real eigenvec- 
tors. We also set 

Vi = (1/11,1/21,1/31)^ (65) 

and 

mi' = iAuAli,,A2iA*2i,,AsiAyf (66) 
for /, /' = 0, 1. The conditions (11) and (13) are then 

A/woi = (67) 



and 

Afwii = Vi . (68) 

Equation (68) has a solution only if Vi is a linear com- 
bination of fii and ^2- Suppose this is indeed the case, 
so that Vi — viifii + V12H2 ■ The equation Af^ = Vi 
has the solutions 

I = (Sii/Ti)/xi + (Si2/T2)At2 + 6M3 , (69) 

where ^3 is an arbitrary real constant. We can set 

wii = (|/lii|M/l2i|MA3in^ = 1 = (6,6,6)^ (70) 

if > for i = 1, 2, 3. Suppose that the ^/s are in fact 
all positive. Then Aji = yj^j exp(i0j), where 0j is real. 
Equation (67) shows that 

Woi = (AioA*i,^20/i;i, ^30/131 )^ 

= C/X3 = C(/i31, ^32, ^33)^, (71) 

where C is an arbitrary nonzero complex constant, and 
so Ajo = C exp{i(f>j) / for j = 1,2,3. Since Ajq 
and Aji are both proportional to exp(i0j), the order pa- 
rameter tpj is proportional to exp(i0j) as well, and we 
may set 0j = for j = 1,2,3 without loss of generality. 
The fij's can be readily obtained from Eq. (14). We con- 
clude that the order parameter for the jth condensate 
component is 

V-, = e/'/'e^^"^* + Cm3j cos(fcix)e-'-i*] (72) 

for j = 1, 2, 3. By changing the zero of time if necessary, 
we can arrange for C to be real and positive. The density 
of the jth condensate component is then 

'^i ^ G^ 7 — ^ cos (kix) 

+2C/i3j cos(fcix) cos(ti;ii), (73) 

which shows that the oscillations of each pair of conden- 
sate components are either in phase or 180° out of phase. 
The phase and amplitude of the oscillations of the com- 
ponents between sublattices 1 and 2 are determined by 
the nature of /^s, the eigenvector of the interaction ma- 
trix AI with eigenvalue zero. 

B. Solutions on a Square Optical Lattice 

For D = 2, the optical lattice is formed by two standing 
light waves with orthogonal wave vectors. We take the x 
axis to lie along ki and the y axis to lie along k2. For 
ki 7^ k2, the optical lattice has rectangular symmetry, 
while for ki = k2 = k, we obtain a square optical lattice. 
In this section, we will study solutions on the square op- 
tical lattice. Solutions on the rectangular optical lattice 
will be discussed in Section IV C. 
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1, One-Component Condensates 

The Ajis must satisfy Eq. (11) for ? = and V = 1, 2. 
They must also satisfy Eq. (12) with / = 1 and V = 2. For 
s = 1, these conditions reduce to ^lo^n = A10AI2 = 
and di{AiiAl2) — 0, respectively. Equation (13) has a 
solution only if pi and gu have the same sign. We assume 
that this is the case. Since \Aii\'^ = Vu/gu > for 
Z = 1, 2, both and Qi must vanish. Wc can arrange 
for All to be real and positive. Ar2 is then imaginary. 
The condensate wave function is thus 

^1 = -, [ei cos(fca:) ± ie2 cos{ky)] e"*"* , (74) 
2 V 511 

where ui = fik^ /2^. This is a stationary solution on the 
square optical lattice with potential 

V = Vi^ -^pi [e\ cos^ [kx) + el cos^ {ky)] (75) 

and timc-indcpcndcnt condensate density ni — —Vi/gn. 
A PC solution therefore exists only if the spatially 
averaged condensate density (ni) is precisely Pi{e\ + 
ei)/(8ffii). 

Although the solution (74) has previously been con- 
structed by Hai et al. [12], its physical interpretation 
has not yet been discussed. The solution is a vortex- 
antivortex lattice [see Fig. 1, Parts (a) and (b)], and so 
the condensate is flowing even though its density is not 
time-dependent. Each square of side A/2 with a potential 
maximum at its center and potential minima at its cor- 
ners is occupied by a vortex or antivortex. The cores of 
the vortices and antivortices are located at the potential 
maxima, where the condensate density is zero. As shown 
in Fig. 1(b), the vortices and antivortices are arranged in 
a checkerboard pattern. 

The mean- field equations of motion (1) are time- 
reversal invariant: if ■0(r, i) is a solution, then so is the 
time-reversed state ■0*(r, —t). Using our method of solu- 
tion, we found both the solution with the upper sign in 
Eq. (74) and its time-reversed version, the solution with 
the lower sign. 



2. Two- Component Condensates 

We have now finished establishing contacts between 
the literature and the solutions obtained using the for- 
malism of Section III. To the best of our knowledge, the 
solutions found from this point on are new. 

The analysis for two-component condensates on a 
square optical lattice runs parallel to that given for two 
components on a one-dimensional optical lattice, and so 
only the final results will be given. If det M 7^ 0, there is 
a solution of the form 



(a) (b) 




FIG. 1: (Color online) Solutions on a square optical lattice: 
(a) Gray scale plot of the square optical lattice potential 
V{x,y). Regions of low (high) potential are shown in black 
(white), (b) The current density for the one-component so- 
lution (74) with the upper sign. The direction (size) of the 
arrows indicates the direction (magnitude) of the current flow. 
This solution is a vortex-anti-vortex array, (c)-(f) Gray scale 
plots of the density of the first component ni{x,y) for the 
two-component solution to the Manakov case with oq — ai, 
p = 1 and 9 = 37r/4. The plots are for times t = T/8 
[Panel (c)], t = 3T/8 [Panel (d)], t = 5T/8 [Panel (e)] and 
t = 7T/8 [Panel (f)]. Regions of high (low) ni are shown in 
black (white). Each plot shows the region with — A < x < A 
and —A < y < X- 



for j = 1, 2, provided that Eq. (13) has a solution for the 
Aji^s. Each of the two condensate components moves in 
a vortex-antivortex lattice in this solution. If the condi- 
tion (49) holds, on the other hand, wc have a solution 
with 



= {\Aji\coakx ± i\Aj2\cosky)e-''^* (76) 



^1 = aoe-'sii""*/' 



(77) 
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and 



V'2 




162 COS ky) 



-i{uj+g2ial/h)t 



(78) 



where ao = ^lo is an arbitrary real constant. For both 

(2) 

Case B and Case C, this is a vahd solution and tp ^ ipi , 
where 



^(2) _ ^-icrgaf,At/h 



ao 



(oi cos kx ± ia2 cos ky)e 



(79) 



There are nonstationary solutions if the problem is fac- 
torizable. For the Manakov Case B, we have an uncount- 

(2) 

ably infinite P-set: the solution if) = P{9)'ipl is valid 
for arbitrary real 6. The densities of the two components 
of the condensate are time-dependent in this solution if 
sm(26) is nonzero. For example, 



ni{x,y,t) 



cos^ 9 + a\ sin^ 6[cos^{kx) + cos^(fcj/)] 
— aoci sin(20)[cos(fca:) cos(a;t) 
-l-pcos(fcy) sin(u;t)], (80) 



where p = ±62/61. The solution with p = —61/62 is 
simply the time-reversed version of the solution with p = 
61/62, and so we may restrict our attention to the case 
p > 0. The external potentials Vi and V2 coincide and are 
equal to V = —jp{el cos^ kx + 63 cos^ ky) . The potential 
minima occur at the points {x,y) = ^{qi,q2)i where qi 
and 92 arc integers and A is the optical wavelength. We 
divide the lattice of potential minima into four square 
sublattices with lattice spacing A, as shown in Fig. 2. 

Although ni and 712 are time-dependent, the total con- 
densate density n does not vary in time, and its maxima 
occur at the potential minima. The time evolution of the 
density of the first component, as described by Eq. (80), 
is illustrated in Fig. l(c)-(f). Let T = 27r/a; be the period 
and suppose for the sake of specificity that oq is positive 
and that 7r/2 < 9 < tt. At time t = T/8, the maxima 
of ni are on sublattice 1 [Fig. 1 (c)] . One quarter period 
later, the maxima of ni are on sublattice 2 [Fig. 1 (d)]. 
They are on sublattice 3 at time t = 5T/8 [Fig. 1 (e)] and 
sublattice 4 at time t ^ 7T/% [Fig. 1 (f)]. Finally, the 
maxima of ni return to sublattice 1 at time t = 9T/8. 
The maxima of 712 also oscillate among the sublattices, 
but the oscillations of n2 lag those of ni by half a period. 

For the FOP Case C, we obtain a P-sct of solutions by 



transforming ■0*^'': explicitly, xp ~ g-«f2At 
is a solution for arbitrary real A. The external potentials 
are Vi = —V2 = V. Suppose for the sake of specificity 
that p is positive. The minima of Vi and the maxima of 
V2 then occur at the lattice of points {x,y) = ^{qi,q2), 
where qi and q2 are integers and A is the optical wave- 
length. We divide this lattice into the same four square 
sublattices as we did for Case B. As time passes, the 



ik — — ik — — ik 



— — ik — <>—tk 



FIG. 2: (Color online) Division of the lattice of potential 
minima into four square sublattices. The lattice of potential 
minima are the vertices of the grid. The vertices in sublattices 
1, 2, 3 and 4 are indicated by circles, squares, diamonds and 
stars, respectively. The origin is at the center of the figure. 



maxima of rii move periodically among the four sublat- 
tices, just as they do for Case B. In Case C, however, the 
oscillations of component 2 are in phase with those of 
component 1, and the total condensate density n varies 
in time as a result. This is analogous to what we found 
for Case C in one dimension. 

For the solutions just discussed, the spatial average of 



the total number density 
Manakov case, while (n) 



i(af -I- al) for the 
cosh(2A)[ag + ^{aj + a|)] 
for the FOP case. In both cases, we obtain solutions 
provided that (n) > ^{af + a^)- 



C. Solutions on Rectangular Optical Lattices 

We now turn to the case in which D ^ 2 and fci 7^/02, 
i.e., to the rectangular optical lattice. No solutions of the 
form of Eq. (8) exist for a single-component condensate 
on a rectangular optical lattice. For a two-component 
condensate, PC solutions are obtained only for the Man- 
akov Case B, and so we will confine our attention to that 
case. Choosing I = and /' = 1,2 in Eq. (22), we ob- 
serve that if Aq is nonzero, Ai and A2 must be parallel. 
Equation (22) with / = 1 and I' = 2 then shows that 
either Ai or A2 must vanish and, hence, Eq. (28) can- 
not be satisfied for both / = 1 and 1 — 2. It follows 
that Aq = 0. We can take the components of Ai to 
be real without loss of generality. Since |Ai| = oi, we 
may set Ai = ai(cos0, sin6')"^. We can arrange for the 
components of A2 to be real through a change in the 
zero of time. Because Ai • A2 = and IA2I = 02, it 
follows that A2 ~ ±a2(— sinfl, cos0)-^. Equation (30) 
shows that = 0, and hence we have the solution given 

by 



i/ji = ai cos 9 cos{kix)e 

=Fa2 sin6'cos(fc2y)6"*'^^* 



(81) 
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and 



iIj2 = aisin9cos{kix)e 

±a2 cos 9 cos(fc22/)e~"^^*, 



(82) 



where the angle 6 is arbitrary. Equations (81) 
and (82) define a P-set of solutions since — 
P{e){ai cos(fcia;)e-*"i*, ±02 cos(fc22/)e-*"^*)^. 

Equations (81) and (82) give a nonstationary solu- 
tion with temporal period T = 2'k/\uj2 — uji\ for < 
9 < tt/2. The time evolution of this solution can 
be understood as follows. The optical potential V = 
— \p[e\ cos^{kix)+e\ cos^(fc2y)] has minima at the points 
(x, y) = 7r(gi/fci, 92/^2), where qi and 52 are integers. Di- 
vide the lattice of potential minima into two sublattices: 
sublattice A with even qi -I- (72 and sublattice B with odd 
91 +92- For the solution given by Eqs. (81) and (82) with 
the lower signs and Q < 9 < 7r/2, the maxima of ni are 
initially on sublattice A, as illustrated in Fig. 3(b) for 
flQ = ai and 9 = -kIA. Half a period later, maxima of 
ni are on sublattice B [see Fig. 3(c)]. The maxima of 
ni are on sublattice A once again at time t = T . The 
second component oscillates between the two sublattices 
in the same way, but its oscillations lag those of the first 
component by half a period. 

In the limit that ki and fc2 coincide, the period of os- 
cillation T tends to infinity and wc obtain a stationary 
solution on the square optical lattice. In this solution, the 
maxima of ni reside on one sublattice and the maxima 
of n2 are on the other. 



D. Solutions on the Simple Cubic Optical Lattice 

For condensates with two or more components on three 
dimensional optical lattices, the set of solutions of the 
form (8) is prohibitively large. Therefore, wc will make 
a number of simplifying assumptions and will limit our- 
selves to giving examples of solutions. The stationary 
solutions we will discuss all have a high degree of sym- 
metry. 

For D = 3, the optical lattice is formed by three stand- 
ing waves with orthogonal wave vectors. The xi axis will 
be taken to lie along k; for Z = 1, 2, 3. We will confine our 
attention to the case in which the three standing waves 
have the same wavelength A, so that the lattice of po- 
tential minima is a simple cubic (SC) lattice with lattice 
spacing A/2. Wc will further simplify the problem by 
restricting our attention to the Manakov Case B and by 
assuming that the ej's coincide. To simplify the notation, 
set fc = fci = ^2 = ^3, w = hk'^ /2m, e = ei = 62 = 63, 
/;(r) = cos(k; • r) for Z = 1, 2, 3, and a = yJp/ge/2. 

From Section III, we know that 



tp = exp(-ig|Aopi/^) 

X [Ao + (Ai/i + A2/2 + A3/3)e-'^"*] 



(83) 



is a solution to the mean-field equations of motion if 
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FIG. 3: (Color online) Two-component condensate on a rect- 
angular optical lattice with Ai = 2A2: (a) Gray scale plot of 
the 2D optical potential V{x,y). Regions of low (high) po- 
tential are shown in black (white), (b)-(c) Density of the first 
component ni{x,y,t) at times t = and T/2, respectively; re- 
gions of high (low) ni are shown in black (white) . (d) Current 
density of the first condensate component at time t = r/4, 
illustrating the flow from sublattice A to sublattice B. The 
direction (size) of the arrows indicate the direction (magni- 
tude) of the current flow. Each of the four plots shows the 
region with — Ai < a; < Ai and — A2 < J/ < A2. 



for I ~ 1,2,3 and 

3fJ(AJ • A2) = 3fi(A* • A3) = 5R(A* • Ai) = 0. (85) 

There is no solution to Eqs. (84) and (85) for s = 1, and 
so we will only consider condensates with two or more 
components. 

For brevity, the lattice of potential minima will be re- 
ferred to as "the lattice." The lattice can be divided into 
eight simple cubic sublattices with lattice spacing A. The 
vector 



f =(/i,/2,/3) 



(86) 



a and An • A; 







(84) 



takes on a different value on each of these sublattices. 
The lattice can also be divided into four body-centered 
cubic (BCC) sublattices. Each of these BCC sublattices 
is the union of two simple cubic sublattices with f 's that 
sum to zero. In Table I, we assign labels to each of the 
eight simple cubic sublattices and to each of the four BCC 
sublattices. These sublattices are illustrated in Fig. 4 and 
will play an important role in our examples. 

All of the stationary solutions we will discuss have at 
least one of the two symmetries we will now define. If, 
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TABLE I: Labeling of the sublattices of the simple cubic lat- 



tice 


f 


SC siiblatticp label 


RCC siiblattire label 




Ui- 


n 


Mil") 
{-1,-1,-1) 


fl — 

yj 





(-1,1,1) 


1+ 


1 


(1,-1,-1) 


1- 


1 


(1,-1,1) 


2+ 


2 


(-1,1,-1) 


2- 


2 


(1,1,-1) 


3+ 


3 


(-1,-1,1) 


3- 


3 



component symmetry. 



1. Two-Component Condensates 

For s ~ 2, all solutions of the form (83) are stationary 
because Eqs. (84) and (85) do not have a solution with 
nonzero Ag. Both examples of solutions we give will 
therefore be stationary solutions. 

One possible solution with four-fold rotational symme- 
try is given by 



^1 



1 

73 



^a(A+/2 + /3)e--*EEV'i™*^ 



(87) 
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FIG. 4: (Color online) Sublattices of the simple cubic lattice: 
The vertices of the grid are sites of the lattice of potential 
minima. The gray sites belong to the SC sublattice 0-1- , while 
the white sites belong to the SC sublattice 0—. The gray and 
white sites together make up the BCC sublattice 0. Sites of 
all eight SC sublattices are labelled at the corners of the cube 
in which x, y and z all lie between and A/2. The sites of 
SC sublattices 1+, 1 — , 2-f , 2—, 3-|- and 3— are colored red, 
cyan, green, magenta, blue and yellow, respectively. 



and 



i'2 



^«(/l+ 6^-^/^/2+6^-/^3)6 



-^ut ^ ^(rot) 



The maxima of the total density are located at the po- 
tential minima. A straightforward analysis reveals that 
ni has its maxima on the BCC sublattice and that the 
maxima of n2 reside on BCC sublattices 1, 2 and 3. This 
solution does not possess component symmetry. 

Component 1 is at rest since the phase of is inde- 
pendent of position. In contrast, the flow of component 
2 is fascinating: it flows in a three-dimensional vortex 
lattice of great beauty, as we will now demonstrate. 

Consider the cube C in which x, y and z range between 
and A/2. Each of the eight corners of the cube belong 
to a different simple cubic sublattice (see Fig. 4). The 
second component of the condensate flows along six of 
the twelve edges of the cube. Speciflcally, component 2 
flows from the site 1+ to the site 3—, and then to sites 
2+, 1—, 3-1- and 2— before returning to site 1+. There 
is no mass current along the remaining six edges of the 
cube. 

The cyclic flow of component 2 suggests that there is 
a vortex line within the cube, and this is fact the case: 
a vortex line has its core along the cube diagonal that 
joins the 0-|- site to the 0— site. To establish this, we will 
begin by considering the behavior of -02 close to the line 
X = y = z. Let 



after a certain lattice translation, the density rij (r) is un- 
changed by a rotation of 90° about the x, y and z axes 
for j = 1, 2, . . . , s, then we say that a solution has four- 
fold rotational symmetry. Note that the lattice transla- 
tion could depend on the condensate component index j 
and could be the null translation. If a solution has four- 
fold rotational symmetry, the x, y and z directions are 
equivalent, and so this symmetry is a type of discrctized 
isotropy. On the other hand, if for each pair (j, j') there 
is a sequence of lattice translations or a series of rota- 
tions of 90° about the x, y oi z axes that maps nj{r) 
onto nji(r), then we say that the solution has compo- 
nent symmetry. Intuitively speaking, the s components 
of the condensate all play the same role in a solution with 



and 



The vectors 
with e'o ~ e- 



1 



1 



2" 2 



^2 = \\^{V-^)-, 



-(f -l-y-Hz). 



(89) 



(90) 



(91) 



'1, 



e' • r, where 



and Gg form an orthonormal triad 
We introduce the new coordinates 
! = 1,2,3. On the line x = y = z, 
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both x'l and x'2 vanish and Xg is arbitrary. Rewriting ip2 
in terms of the new coordinates and expanding for small 



2. Three- Component Condensates 



and x'r,, we obtain 



^"2 



-ka sin 



(92) 



Equation (92) shows that there is a vortex line with its 
core along the line x ^ y ~ z. The direction of the vector 
Gg and the direction of the current flow around the vortex 
core arc related by the right hand rule; for brevity, we 
will say that the vortex line is oriented along the vector 

We have just shown that there is a vortex line within 
the cube C with its core along the cube diagonal, and 
that the vortex line is oriented along the vector 63. To 
determine the nature of the flow throughout space, first 
note that '02 is invariant under the transformation x 
—X, i.e., it is invariant under reflection about the y — z 
plane. Naturally, -02 is also invariant under the reflections 
y —y and z —z. These invariances give the flow 
within the cubical region in which x, y and z range from 
— A/2 to +A/2. Since tp2 is invariant under the lattice 
translation r r + (qix + q2y + q3z)X for all integers qi, 
(j2 and ^3, the nature of the flow in the whole of space 
can now be inferred. 

The following picture emerges from this analysis. 
There is a vortex core along every line in the BCC sub- 
lattice that joins an infinite chain of nearest-neighbor 
sites. At these sites, the density of the second compo- 
nent is at a maximum, its current density is zero and 
the potential is at a minimum. Each vortex core also 
passes thorough a chain of neighboring potential max- 
ima which alternate with the potential minima. At the 
potential maxima, the density of the second component 
of the condensate 712 is zero. In this way, the energy of 
the vortex array is minimized. Every vortex line is ori- 
ented along one of the following four vectors: x + y + z, 
X — y — z, —X + y — z or —x — y + z. 

The solution with 



/2 + e 



-i-iT/4 



(93) 



(rot) 



For s = 3, the stationary solution given by tfji — 

^2 = V0V'2''°*^ and 03 = Vl~~0V'2''°*^ with < < 1 
has four-fold rotational symmetry but does not have com- 
ponent symmetry. Next, consider the stationary solution 
with 



- \a U + ./2 + h \fj 



(95) 



for i = 1,2,3. The maxima of Uj are on the jth sublattice 
and each component of the condensate is at rest. This so- 
lution does not have four-fold rotational symmetry since 
the Xj direction is special for condensate component j. 
However, it does have component symmetry. To see this, 
consider an arbitrary pair of indices {I, V) with I ^ I'. 
Let I" be the integer belonging to the set {1,2,3} that 
differs from both I and I' . A 90° rotation about the xi" 
axis interchanges // and //', and so maps n; onto n//. 

A nonstationary solution that illustrates just how com- 
plex the solutions for s = 3 can be is given by 



01 



V2 



(/i+/2 + */3)e-'("+'')*, 



02 



and 



0'3 



V2+(-/i + /2 + i/3)e- 



V2 - (-/i + /2 + if3)e- 



(96) 



(97) 



(98) 



and 



where fl = ga^/h. The density of the first conden- 
sate component is time-independent and its maxima are 
on BCC sublattices and 3. The density maxima of 
component 2 are on simple cubic sublattice 1+ ai time 
t = ^ tan~^(l/2) = r, on simple cubic sublattice 2+ at 
t = r/2 — T, on simple cubic sublattice 1— at t = T/2 + t, 
and are on simple cubic sublattice 2— at time t = T ~ t. 
At time T + t, the maxima of n2 have returned to simple 
cubic sublattice 1+. The motion of condensate compo- 
nent 3 is identical to that of the second component, ex- 
cept that the oscillations of ns lag those of 7^2 by half a 
period. 



02 



V2 



(/i 



-i7r/4 r 



(94) 



has component symmetry but not four-fold rotational 
symmetry: The x direction is not equivalent to the y and 
z directions, and ni and n2 differ only by a translation 
through the distance A/2 along the x axis. The maxima 
of 111 are on BCC sublattice 0, while the maxima of 712 
are on BCC sublattice 1. 

A natural question to ask is whether, for a given s, 
there is a solution with both four-fold rotational and 
component symmetries. The answer to this question is 
"no" for both s — 2 and 3, as we show in Appendix B. 



3. Four- Component Condensates 

The solution space for four-component condensates is 
very large. To see this, consider an arbitrary set of or- 
thonormal vectors with real components in four dimen- 
sions, {cq, ci, £2, £3}. Setting Aq ~ aoco and A; ~ aii 
for I = 1,2, 3, we obtain a solution to Eqs. (84) and (85) 
for arbitrary nonnegative real numbers oq. One such so- 
lution is given by 



1 



[-ao + a(/i + h + /3)e-'"*] e"''^'^'*/'' (99) 
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and 



i'j = 2 ^ "'^■^^ + f2 + h~ 2/j-i)e" 

y-p-igo-lt/n 



(100) 



for j = 2, 3, and 4. 

For two- and three-component condensates, no station- 
ary solution of the form (83) has both four-fold rotational 
and component symmetries. Such a solution does exist 
for four-component condensates, however. The solution 
is given by Eqs. (99) and (100) with oq = 0. In this solu- 
tion, the maxima of rij are on sublattice j — 1 and each 
of the four condensate components is at rest. To see that 
the solution has four-fold rotational symmetry, note that 
ni is unchanged by a rotation of 90° about the x, y and 
z axes. After a translation through ^Axj-i, the density 
Uj becomes rii and so is unchanged by a rotation of 90° 
about the x, y and z axes for j = 2, 3, 4. 

As we have seen, Uj can be mapped onto ni by a prim- 
itive lattice translation for j = 2, 3 and 4. On the other 
hand, ni is mapped onto rij' by a translation through 
^Xxj'-i for j' = 2, 3, and 4. It follows that, for each 
pair there is a sequence of at most two primitive 

lattice translations that carries rij onto riji , and so the 
solution has component symmetry. 

For ao > 0, Eqs. (99) and (100) describe a nonstation- 
ary solution. In this solution, the density maxima of the 
first component of the condensate are on the simple cubic 
sublattice 0— at time < = and are on the simple cubic 
sublattice 0-t- att^ T/2. For j = 2, 3 and 4, the maxima 
of rij are on the simple cubic sublattice j+ at time t ~ 
and are on the simple cubic sublattice at t ~ T/2. 
The entire condensate returns to its initial state at time 
t = T. 



The elements gjji of the interaction matrix must also 
be put into dimensionless form, but first we note that 
they may be renormalized, depending on the dimen- 
sionality of the optical lattice. If there is narrow har- 
monic transverse confinement for one-dimensional and 
two-dimensional optical lattices, the appropriate forms 
in all dimensions are [40] 



(3) 



(2) 



\ m 



1/2 



a 



n' ' 



(104) 



and g^^^i = 2huj±ajj' , 



where ajji is the low-energy s-wavc scattering length for 
species j and j' , the superscript on gjj' is D, the dimen- 
sionality of the optical lattice, and the confining poten- 
tials are characterized by the angular frequencies coz in 
two dimensions and uj±_ in one dimension. 

We may choose the elements gjj' of the dimensionless 
form of the interaction matrix to be typically of order 
one, so that the elements of the scattering-length matrix 
and the dimensional interaction matrix decompose as 



' = agjj> and 5 / ^ g^ >g,y 



(105) 



where a and g*-^-* are scalar, dimensional factors. The 
latter is the same as the g appearing in Eq. (19), but 
with the possible need for renormalization in lower di- 
mensions explicitly indicated by the superscript. The 
dimensionless form of is 



(D) 



Eox§ ■ 

Its values, which follow from Eqs. (103)"(106), are 



(106) 



V. NONLINEAR STABILITY 
A. Dimensionless mean-field equations 

Our numerical investigations of the stability of selected 
solutions to the mean-field equations (1) are performed 
with a dimensionless form of the equations, using dimen- 
sionless position, time, and potential-energy variables, 

^ = r/.To, T = t/to, and Vj=Vj/Ea, (101) 

defined in terms of units xq, toj and Eq. For optical 
potentials of the form (3), we also define dimensionless 
potential-strength coefficients 



V,i^V,i/Eq. 



(102) 



The length and time units, xo and are related to the 
energy unit Eq by 



Xq 



h/^/mEQ and = h/Eo . (103) 



;;(3) 
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47r 



ymEoa , 



inmuj. 



1/2 



and g^ ' =2. -—uj_i_a . 

V En 



(107) 

We absorb the square root of this dimensionless scale 
factor into the order parameter, which then takes the 
dimensionless form 



(108) 



Thus, for a PC solution, the coefficients appearing in ^pj 
are related to those in Eq. (8) by 



A 



(109) 



The normalization of the dimensionless order parameter 
is given by 



E J\Mtr)\'d^^N~g(''^ Vr, 



(110) 
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wherein we see that g'^^^ plays a role equivalent to the 
number of particles, N. For a PC solution, the dimen- 
sionless mean number density is then 



(n) 



X? (n) 



~i)E(i^op4Eii 



(111) 



where (n) is given in Eq. (17). 

Finally, the mean-field equations (1) take the dimcn- 
sionless form 



' dr 



(112) 



B. Details of the calculations 

Our numerical stability tests use Eq. (112), propagat- 
ing a specified initial condition ■0^(^,0) forward in time 
via a fifth-order Runge-Kutta algorithm with adaptive 
step-size control [41], the spatial derivatives being calcu- 
lated in wave- vector space via a pseudo-spectral method. 

We perturb the solution by adding some white noise to 
the initial condition before beginning the time propaga- 
tion. This is accomplished for each component of the or- 
der parameter by adding to the real and imaginary parts 
of each of its Fourier components a random number from 
a uniform distribution in the range ±0.5 x 10~^ times the 
modulus of the largest Fourier component. 

We work within a spatial cell comprising four periods 
of the optical lattice (two optical wavelengths) in each 
of the D dimensions, applying periodic boundary condi- 
tions to that cell. The spatial grids contain rigiid = 128 
points for one-dimensional cases and rigi-id = 32 points in 
each dimension for two-dimensional cases. Thereby we 
are able to test the stability of solutions against pertur- 
bations having wavelengths ranging from 2A/ngiid to 2 A, 
where A is the optical wavelength. 

As a measure of the instability of a component of a 
solution at a particular instant of time, we use the vari- 
ance of its Fourier power spectrum relative to that at 
T = [39, 43], 



\ 



2^[jK0)f 



(113) 



Here fj{K,T) 
form of ipj (^ 



|2 7 



= |0j(K, t)| , 4>j{K, t) is the Fourier trans- 
r), and Ki = 

A solution is deemed to have reached the onset of in- 



stability when each of the Uj has exceeded 0.1 at least 
once. We find that this criterion correlates nicely with 
the visual onset of instability in the graph of the density 
and works well for a range of solution types and potential 
strengths. 



The time of onset of instability can be sensitive to 
many details, including the amount of added noise, the 
resolution of the spatial grid on which the solution is 
represented, and even details of the generation of the 
random deviates and the algorithm used to perform the 
fast Fourier transforms, particularly when the solution is 
stable for long times. As well, we expect the lifetime of 
an experimentally produced condensate to vary with the 
level of noise present. Consequently, one should not in- 
fer from our graphs of instability-onset time vs. solution 
parameters that the times represent literal lifetimes that 
would be observed in any particular experiment. 

However, as will become clear below, there is a fairly 
well-defined boundary between unstable solutions and 
stable solutions, beyond which the instability-onset times 
increase extremely rapidly. The locations of those bound- 
aries are largely insensitive to details of the calculations. 
We therefore expect the parameter boundaries delimiting 
numerically stable solutions to be experimentally mean- 
ingful, in the sense that within the stable regions ob- 
served lifetimes should be at least of order one second. 



C. Results of the calculations 

To make our results more concrete, we have chosen 
typical values for the laser wavelengths, A = 800 nm in 
all directions, the trap frequencies, = 27r x 100 Hz 
and = 27r x 200 Hz, the atomic mass, m = 87 u, and 
the s-wave scattering length, a = 55 A. We will refer to 
these below as the "system parameters." 

The results are displayed primarily in recoil units, set- 
ting 



2l2 



Eq — Ej- 



2 m 



0.172 /iKxfcs (114) 



and 



to = tR = Ti/Er « 4.44 X 10~^ s , 



(115) 



where fcz, is the laser wave number, and Ub is Boltz- 
mann's constant. The numerical values are obtained from 
our chosen system parameters above. 

We present below selected numerical stability analyses 
for one to three condensate components in one and two 
dimensions. Three-dimensional cases are not included, 
for they are too computationally demanding at this time. 

It is straightforward to transform the system-specific 
values shown in the figures below, the potential strengths, 
the instability-onset times, and the numbers of particles 
per well, to values appropriate for alternative choices of 
the system parameters. We will elaborate on this point 
in Sec. VD, following the presentation of the results. 
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FIG. 5: (Color online) Instability-onset times ti for a one- 
component PC solution on a square optical lattice. The 
dimensionless potential-strength parameter V/Er and the 
corresponding particle density are shown on the horizontal 
axes. The solution is stable at both low and high potential 
strengths. 



1. One component on a square lattice 

We choose the dimensionless interaction parameter 
and the dimensionless potential coefficients Vji to be 

.911 = 1.0, Vh^V/Er, and Vvi^V/Er, (116) 

where the potential-strength parameter V can be var- 
ied. Then the dimensionless solution corresponding to 
Eq. (74) with the upper sign has coefficients 



0, ^11 



Vii , and Ai2 



(117) 



Because the spatially constant term is required to vanish, 
the particle density, Eq. (17), is uniquely determined by 
the potential-strength parameter. 

The instability-onset times for this solution are shown 
as a function of V in Fig. 5. Two time scales are in- 
cluded, showing both recoil times and milliseconds, with 
a maximum propagation time of several seconds. The 
top scale shows the particle density in particles per well 
corresponding to the potential-strength parameter shown 
on the bottom scale. 

While the solution becomes unstable in just a short 
time over much of the range of potential strengths shown, 
it is stable for sufficiently weak potentials. Much more 
surprisingly, it is also stable for sufficiently strong po- 
tentials. To test whether the solution becomes unstable 
again for potentials stronger than that at the boundary 
near 3.5i?_R, we performed propagations to t « 1542^^ « 
68 ms of solutions having V /En as high as 48, finding no 
recurrence of instability. This extends well into the Mott 
insulating regime, beyond the point of physical relevance 
of the mean-field equations, as discussed in Sec. I. 
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FIG. 6: (Color online) Regions of stability and instability for 
a PC solution having two components on a one-dimensional 
optical lattice. The points denote calculated boundaries be- 
tween stable and unstable behavior, the lines connecting them 
serving as guides to the eye. The inaccessible region has no 
PC solutions of the form (119). 



2. Two components on a one- dimensional lattice 

Here we study the stability for the Manakov case, in 
which the dimensionless interaction matrix has rank one 
and has all elements equal to one. The potential coeffi- 
cients we choose are 



(118) 



and the coefficients of the dimensionless solution corre- 
sponding to Eq. (54) are 



AiQ = a , A20 = , 

All = , and A21 = yVii, 



(119) 



where a is a free parameter. The components of this 
solution are then mixed using P{6) of Eq. (56) with 9 = 
7r/4 to produce a nonstationary solution. 

Since the spatially constant term in the resulting solu- 
tion is not required to vanish, we can vary the parameter 
a to control the particle density independently of the 
potential-strength parameter, giving a two-dimensional 
domain in which to investigate the stability of the so- 
lution. As is clear from Fig. 5, the boundaries of the 
stable regions are approximated well by the positions at 
which the instability-onset times have exceeded about 
three hundred times tji . Consequently, we have used that 
as a threshold to define those boundaries for the present 
case in scans over the potential-strength parameter at 
several fixed values of the particle density. The resolu- 
tion of the potential grid was O.WEr, easily adequate for 
the graphical delimitation of the stable regions. The re- 
sults are shown as a map of stable and unstable regions 
in Fig. 6. 
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The boundary of the region marked "inaccessible" cor- 
responds to the vanishing of the coefficient Aiq of the con- 
stant term in the solution. Within that region it is impos- 
sible to construct a PC solution of the chosen form, since 
Eq. (34) cannot be satisfied. The region marked "stable" 
is that portion of the parameter space where the criterion 
for stability is satisfied, and the region marked "unstable" 
corresponds to parameters for which the instability-onset 
time falls below the threshold. As in the two-dimensional 
case shown in Fig. 5, the solution becomes stable in the 
weak-potential limit. However, in striking contrast to the 
two-dimensional case, there is no evidence of a second re- 
gion of stability at high potential strengths. 

In order to verify that apparent absence of stabil- 
ity for deep potentials, we performed calculations along 
the boundary of the inaccessible region with V/En as 
high as 48, well beyond the highest shown in Fig. 6, 
finding only monotonically decreasing instability-onset 
times [44]. This confirms the observation in Fig. 6 that 
there is no additional stable region at high potential 
strength. 



3. Two components on a square lattice 
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FIG. 7: (Color online) Instability-onset times for a two- 
component PC solution on a square optical lattice. The 
dimensionless potential-strength parameter V/Eh and corre- 
sponding particle density are shown on the horizontal axes. 
The solution is stable in three regions, having low, high, and 
intermediate potential strengths. 



As we did in one dimension, here we study the Man- 
akov case, with all elements of the dimensionless interac- 
tion matrix equal to one. Now there are four potential 
coefficients, which we choose to be equal: 



1^11 = 1^12 = V21 = V22 = V/Er . 



(120) 



The coefficients of the solution corresponding to Eq. (79) 
are 



A 



10 



a. 



A 



20 



0, 



^11=^12=0, A21 = v^i 



A. 



and 



(121) 



where a allows us to set the particle density. Once again 
we mix the components using P{0) with 9 = 7r/4 to ob- 
tain a nonstationary solution. 

Instability-onset times for this solution with a fixed at 
one and varying potential strengths are shown in Fig. 7, 
where the conventions are similar to those used for the 
single-component case in Fig. 5. As in that case, re- 
gions of stability occur at both low and high potential 
strengths, but now a rather striking additional region of 
stability appears at intermediate potential strengths, just 
below 10 recoil energies. 

To explore this behavior in more detail, we map out 
regions of stability in the plane of particle density and 
potential strength in Fig. 8 using the same strategy ap- 
plied in the one-dimensional case in Fig. 6. Three areas 
of stability are clearly evident, though the central one is 
somewhat narrower than the others. The fine black line 
running parallel to the boundary of the inaccessible re- 
gion is the track in the parameter-space plane followed 



by the graph shown in Fig. 7. We extended the search for 
renewed instability along this line to V/Er = 48, limit- 
ing the propagation time to t 1542tR ~ 68 ms, finding 
no evidence of further instability beyond the crossover 
into the stable region near V/Er = 16. 
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FIG. 8: (Color online) Regions of stability and instability for 
the PC solution Eq. (121), having two components on a square 
optical lattice. The points denote calculated boundaries be- 
tween stable and unstable behavior, the lines connecting them 
serving as guides to the eye. The inaccessible region has no 
PC solutions of the form Eq. (121). The fine line parallel to 
the boundary of the inaccessible region is the track followed 
by the graph of instability-onset times in Fig. 7. 
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V/Er 

FIG. 9: (Color online) Regions of stability and stability for 
a PC solution having three components on a one-dimensional 
optical lattice. The points denote calculated boundaries be- 
tween stable and unstable behavior, the lines connecting them 
serving as guides to the eye. The inaccessible region has no 
PC solutions of the form Eq. (124). 



4. Three components on a one- dimensional lattice 

We have also tested the stability of a class of PC so- 
lutions having three components. The dimensionless in- 
teraction matrix then has nine elements, all ones in the 
Manakov case, and rank equal to one. The potential co- 
efficients are all chosen to be the same: 

Vii = V21 = %i = V/Er . (122) 

The coefficients of the components of the dimensionless 
solution are all of equal magnitude, those of the spatially 
constant term being real: 

Aio = A20 = iao = a , (123) 

with those of the space-dependent part chosen to have 
the phases of the cube roots of unity: 




and 

(124) 



D. Alternative choices of system parameters 

Those aspects of the results presented above that are 
dependent on the system parameters chosen in Sec. V C 
are readily transformed to values corresponding to al- 
ternative choices of those parameters. Obviously, the 
potential-strength parameter V is trivially obtained by 
multiplying the dimensionless value V/Er by the recoil 
energy corresponding to any desired set of system param- 
eters, and the instability-onset time ti is easily converted 
by multiplying it by the ratio t'^^/tR^ of the time units 
t'^ corresponding to the alternative parameters and Ir 
corresponding to our chosen system parameters. 

The average number of particles per well is just the 
mean density times the well volume, 

^weii- f 2 ) • 

From the dimensionless density given in Eq. (Ill), we 
see that this can be expressed in terms of dimensionless 
parameters as 

3 — 1 ^ — 1 

(126) 

with the length unit xq set to the recoil length xr = 
Ti/^mER. 

For all of our stability figures, the dimensionless 
potential-strength parameter V/Er is a free parameter, 
and it determines the values of the dimensionless coeffi- 
cients Aji having I > 0. For Figures 6, 8, and 9, there is 
one additional free parameter, a, and it determines the 
values of one or more of the coefficients Aj^. Thus, for 
any given abscissa in Figure 5 or 7, or abscissa and ordi- 
nate in Figure 6, 8, or 9, the sum in Eq. (126) is fixed, 
and the value of iVwcii corresponding to an alternative 
choice of system parameters can be obtained from that 
shown on the graph by simply rescaling the prefactors: 

where the unprimed quantities correspond to our choice 
of system parameters, and the primed quantities to some 
alternative choice. 



VI. CONCLUSIONS 



The map of stable regions in the space of the free pa- 
rameters of this solution is shown in Fig. 9, which is vi- 
sually indistinguishable from its two-component analog, 
Fig. 6, having a region of stability where the potential is 
sufficiently weak. In fact, the coordinates of all the points 
plotted on the graph arc identical, within the resolution 
of the scans. 



In this paper, we made a comprehensive study of 
potential-canceling (PC) solutions for an s component 
Bose-Einstein condensate in in a D-dimensional optical 
lattice. Studies of specific cases with small s and D, es- 
pecially in one spatial dimension, have appeared in the 
literature. Our work brings these previous studies to- 
gether, generalizes to arbitrary s and _D, and provides 
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intriguing new solution types and novel physical inter- 
pretations. 

Currently, there is a great deal of interest in the 
a Berezinskii-Kosterlitz-Thouless phase in Bose-Einstein 
condensates at intermediate temperatures in 2D [45] . In 
such a phase, vortex-anti-vortex pairs become bound to- 
gether, in contrast to the free vortex proliferation which 
occurs at high temperatures. However, this phase is re- 
stricted to a truly 2D system, which is difficult to achieve 
experimentally. We have shown that an optical lat- 
tice stabilizes vortex- anti- vortex pairs in the quasi-two- 
dimensional case, and that the lattice causes the array to 
be tightly packed. 

Not only have we presented multicomponent gener- 
alizations of 2D vortex-anti-vortex arrays, but we have 
also generalized them to three dimensions. In 3D, we 
constructed a solution in which one condensate compo- 
nent forms a lovely and complex three-dimensional ar- 
ray of intersecting vortex lines. As a part of our study 
of PC solutions in 3D, we gave a thorough treatment of 
the most highly symmetric solutions for condensates with 
two, three and four components. Our formalism can also 
be used to gain insight into complex systems now exper- 
imentally available, such as five-component condensates 
in 3D optical lattices. 

We studied the stability of PC solutions numerically 
in ID and 2D for one-, two and three-component con- 
densates. We found three main results: (1) potential- 
canceling solutions tend to become stable as the poten- 
tial strength is reduced; (2) there is a remarkable differ- 
ence between the one-dimensional and two-dimensional 
solutions, in that the latter are also stable for deep po- 
tentials; and (3) for two-components in a square optical 
lattice, there is a fascinating third region of stability for 
intermediate-strength potentials. We found no evidence 
of stabilization at high potential strength in one dimen- 
sion. 

Finally, we mention that the possibility of experimen- 
tally realizing vortex-anti- vortex arrays in a 2D lattice is 
provided for in the recent experiments of Sebby-Strabley 
et al. [46]. In those experiments, two polarizations of 
the lasers used to create the optical lattice potential are 
manipulated to create lattices that can be dynamically 
controlled on a site-by-site basis. In this way, one can 
imagine creating an array of small "propellers" to stir 
up vortex-anti-vortex pairs. The parameter ranges in 
which such a procedure would lead to stable structures 
were determined in our numerical studies. To manipu- 
late multicomponent condensates, one can imagine more 
advanced versions of such an experiment, in which the 
fact that different hyperfine components "feel" different 
lattice strengths for a given optical wavelength A can be 
used to one's advantage. 

We thank B. Deconinck, J. N. Kutz, and J. N. Roberts 
for useful discussions. LDC's work was supported by the 
National Science Foundation under Grant PIIY-0547845 
as part of the NSF CAREER program. 



APPENDIX A: PROOF THAT THE A/S CAN BE 
RESCALED TO HAVE UNIT MODULUS 

For Special Case A, the equations of motion are fac- 
torizable and are given by Eq. (35). The elements of the 
interaction matrix are gjji = agXjXji and the optical po- 
tentials are Vj = XjV. Consider an associated "normal- 
ized" problem that is also factorizablc. In this normal- 
ized problem, the elements of the interaction matrix are 
gjji = agXjXji and the optical potentials are Vj = XjV, 
where Xj = Aj/|Aj| has unit modulus. Suppose we have 
a solution 

s 

to the normalized problem. We can then construct a cor- 
responding solution to the original, unnormalized prob- 
lem as follows. We let Aji = Aji/y/\Xj \ and Qj = \Xj\i}j, 
and define i/> through Eq. (8). if) is then a solution to 
Eq. (35). Moreover, \ipj\'^ = [i/ijlVl^il- We see that 
for each solution ^p of the normalized problem, there is 
a corresponding solution ^p to the original, unnormalized 
problem, and that the density of the jth condensate com- 
ponent simply differs by the constant factor |Aj|^^ in the 
two problems. As a result, we may assume without loss 
of generality that the A^ 's all have unit modulus. 

The length of Aq is a free parameter at this point. 
However, if the average total density (n) is specified in 
the original, unnormalized problem, then Eq. (17) gives 

J=l ^ \ 1=1 I 

If Eq. (A2) has a solution, it fixes the value of |Ao|^. We 
conclude that |Ao|^ is determined if (n) is given. 

APPENDIX B: PROOF THAT THERE ARE NO 
SOLUTIONS WITH BOTH FOUR-FOLD 
ROTATIONAL AND COMPONENT 
SYMMETRIES FOR TWO- AND 
THREE-COMPONENT CONDENSATES IN 3D 

It was stated in Sec. IV D that there are no solutions 
with both four-fold rotational and component symme- 
tries in 3D if the condensate has two or three compo- 
nents. Our proof is as follows. Let s be 2 or 3. The 
condensate order parameters are 

i>,= (^A,,]^e-^-\ (Bl) 

where j ranges from 1 to s and /; = cos(k; • r). The 
densities are 

3 

".=El^^'l'/' + 2 E 'St{AflA*y)h!v (B2) 

1=1 l<i<i'<3 
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If rij is to possess four-fold rotational symmetry, the co- 
efficients of /f , /I and /| must be the same. Thus, \Aji\ 
must be independent of L If the solution is to have the 
component symmetry, on the other hand, \Aji\ cannot 
depend on j. It follows that l^j/p = a^/s for all j and I. 

By choosing a phase, we can arrange for Aji to be real 
for j e {l,2,...,s}. Let Aj2 = ae'^j/yi and Aj^ = 
ae^^' I y/s, where ctj and /3j are real. Then 

n, = ^{A2 + /| + /| + 2[cos(a,)/i/2 

+ cos(«, - i3j)hh + cos(/3j)/3/i]} . (B3) 

For rij to have four-fold rotational symmetry, we must 
have 

I cos I = I cos(Q!j — /?j ) I ~ I cos /?j | , (B4) 

while the condition 

I cos q;i I = . . . = I cos cts I (B5) 

must be satisfied if the solution is to have component 
symmetry. 



Equation (B5) must be reconciled with the condition 
5R(At • A2) = 0, i.e., 

S 

^ cos a.j ^ 0. (B6) 

This is not possible for s = 3, and so there is no solution 
with both four-fold rotational and component symme- 
tries in that case. 

The case s = 2 requires further analysis. Equa- 
tion (B6) gives a2 = TT + aiai, where cti = ±1. Similarly, 
the condition 3?(A3 • Ai) = implies that /32 = tt + (T2/3i, 
where (T2 = ±1. If cti = (T2, the condition ■ A3) = 

becomes cos(q;i — f3i) = 0. Referring to Eq. (B3), we 
observe that if ni is to have four-fold rotational symme- 
try, cos ai and cos (3i must vanish as well. This is not 
possible. If (7i = — 1T2, on the other hand, the condition 
3fi(A2-A3) = becomes cosai cos/3i = 0. Equation (B3) 
shows that if ni is to have four-fold rotational symmetry, 
it is required that cosai = cos/3i = cos(ai — = 0, 
which is an impossibility. Wc conclude that there is no 
solution with both four-fold rotational and component 
symmetries for two components. 
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